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In this paper we present a new model for modeling the diffusion and relative dispersion of particles 
in homogeneous isotropic turbulence. We use an Heisenberg-like Hamiltonian to incorporate spatial 
correlations between fluid particles, which are modeled by stochastic processes correlated in time. 
We are able to reproduce the ballistic regime in the mean squared displacement of single particles 
and the transition to a normal diffusion regime for long times. For the dispersion of particle pairs 
we find a t 2 -dependence of the mean squared separation at short times and a i-dependence for long 
ones. For intermediate times indications for a Richardson t 3 law are observed in certain situations. 
Finally the influence of inertia of real particles on the dispersion is investigated. 
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I. INTRODUCTION 

Relative dispersion of particles in turbulent flows is of 
key importance in a variety of natural and industrial pro- 
cesses, ranging from the spreading of clouds in the atmo- 
sphere or pollutants in the ocean to mixing in pneumatic 
conveyors or production of nanoparticles in flames. There- 
fore the diffusion of particles and the relative dispersion 
of pairs of particles in homogeneous isotropic turbulence 
have been investigated intensively p~r[7] . Although consid- 
erable progress in understanding these processes has been 
made [51 [7J , we still lack the fundamental understanding 
of the observed phenomena. 

Recent advances in direct numerical simulations (DNS) 
of isotropic turbulence [8] made it possible to investigate 
turbulent dispersion [5] at increasingly high Reynolds 
numbers. But because of their extremely high demand of 
computer resources, DNS calculations are unfortunately 
not applicable in most situations. To circumvent this dif- 
ficulty many different techniques for turbulence modeling 
[5] have been invented. 

In this paper we present a new model, which is especially 
useful for the investigation of particle dispersion at high 
Reynolds numbers. Based on a system of stochastic dif- 
ferential equations (SDEs) introduced by A. M. Reynolds 
[TOUT!?] that describes the fluid by a set of tracer parti- 
cles, we use a Heisenberg-like Hamiltonian to incorporate 
spatial correlations between fluid particles. This enables 
us to generate a three dimensional turbulent velocity field 
for which the spatial structure functions show the desired 
behaviour. We then investigate the dispersion of single 
fluid particles by measuring the mean squared displace- 
ment (MSD) as well as the relative dispersion of pairs 



of fluid particles through their mean squared separation 
(MSS). Finally we replace the fluid particles by real heavy 
particles and investigate the influence of the inertia of 
these particles on the MSD as well as the MSS. 

Turbulent dispersion was first investigated by Tay- 
lor [T3], who studied the diffusion of single particles in 
isotropic turbulence. Later Richardson pQ, Obukhov [2J, 
and Batchelor [3] extended the investigation to multiple 
particles. Many measurements of the spread of clouds 
or puffs, as well as the separation of rather large tracer 
particles have been performed (see e.g. Ref. [14] . pp. 556- 
567). Recent advances in experimental techniques made 
it possible to track an increasing number of small trac- 
ers at high Reynolds numbers in laboratory flows (see 
e.g. Ref. [7]). In three dimensions experiments are of- 
ten performed in a tank filled with water and turbulence 
is generated e.g. by counter-rotating baffled disks [15] 
or rotating propellers located in the corners of the tank 
[16] . In two dimensions experiments have been performed 
with thin layers of conducting fluids where turbulence is 
generated through electromagnetic forces [XT] . Simula- 
tions most often use DNS methods in three dimensions, 
where the Navier-Stokes equations are solved directly and 
turbulence is generated by a spectral forcing technique 
[SI [TS] . Kinematic simulations [TH] are another technique 
frequently used in simulations. Most of these publications 
focus on the investigation of the MSS of particles as well 
as the pair-separation probability density functions. More 
detailed introductions can be found e.g. in the reviews of 
Sawford [5] and Salazar & Collins [7J. 

II. SPATIAL CORRELATIONS BETWEEN 
FLUID PARTICLES 
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The spatial and temporal structure of homogeneous and 
isotropic turbulence is highly non-trivial. Of particular 
interest are the spatial velocity structure functions. Let 
u (x, t) be the fluid velocity at time t and position x. 
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The (longitudinal) velocity structure functions are then 
defined as 



(1) 



D n (r) = ( (u(x + r,t)-u(x,t))-r 



where r is a vector connecting two points in the fluid 
field that are separated at a distance r = |r|. According 
to the Kolmogorov similarity theory [14 these structure 
functions depend only on the mean energy dissipation 
rate (e) and the distance r, i.e. 



D n (r) = C n ({s) ■ r) 



<n 



(2) 



Further by Kolmogorov's second hypothesis the exponents 
Cn should be given by 



However, several measurements [2U] show that for n > 3 
the exponents are smaller than the Kolmogorov values. 
Most commonly studied [2T] is the second order structure 
function, for which it can be shown [5] that 
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r <C rj 



D 2 (r)= {C 2 ((e)-r) 2/3 V « r « L , (4) 
2cr 2 L«r 

where v is the kinematic fluid viscosity, a\ the variance of 
the fluid velocity, r\ = [v 3 / (e)) 1 ^ the Kolmogorov length 
scale and L = ct„/ (e) the turbulence length scale. The 
value of the Kolmogorov constant [5] is Ci = 2.0. 

To generate a turbulent velocity field we start with a 
model introduced by Reynolds |10j . This model describes 
well the measured [22] velocity and acceleration distri- 
butions of Lagrangian tracer particles in fully developed 
turbulence by modeling the motion of non-interacting 
fluid tracer particles that are self-correlated in time. It 
consists of a set of stochastic differential equations (SDEs) 
for the logarithm of the local energy dissipation rate 
X = In (ej (e)) as well as equations for one component of 
the acceleration, velocity and position of a fluid particle: 

d X = ~ (X ~ (x^T-'dt + fagX 1 ^ (5a) 
T7 1 + t- 1 - a:i dC ^] a t dt - T^t-^dt 



dot 



Q *l e dt 



2al {T^+t^TzH^db 

dut = atdt 
dx t = Utdt. 



(5b) 

(5c) 
(5d) 



The variance of x m Eq. ( 5a ) is approximated [23] by 



cr~ = —0.354 + 0.289 log R\, where the Reynolds number 

at the Taylor microscale is given by R\ = WlBa^L/ (e). 
The mean value is given by (x) = — 0.5ct^, and the relax- 
ation time scale by T x = 2g\ /(Co (e)). Further in Eq. (5b) 



the energy-containing time scale (sometime also called in- 
tegral time scale) is given by TL = 1<j\I(CqE), the energy- 
dissipation time scale by t v — C ^ 1 / 2 /(2a e 1 / 2 ), and the 
conditional acceleration variance is a 1 , = an£ 3/ ' 2 z/ 1 / 2 . 

at|e u 

The two universal Lagrangian velocity structure constants 
are given as ao = 3.3 and Cq — 7.0. Finally d£i and d^2 
are two independent Wiener processes, i.e. Gaussian dis- 
tributed random numbers with zero mean and variance 
dt. 

In Ref. [23] we already used this model to investigate 
the mixing of heavy particles in a turbulent channel flow 
due to intrinsic fluctuations in the fluid velocity. There 
we "glued" one fluid tracer particle to every real particle 
and used the resulting velocity vector u t in an empirical 
drag law to get a force acting on the real particle. The 
equations of motion of the fluid particles were integrated 

and their positions were deter- 



(3) according to Eqs. (5a )-( 5c 



mined by the positions of the real particles. The biggest 
deficiency of this model is that there are no spatial correla- 
tions between the fluid particles. Calculating the velocity 
structure functions ([T]) in this model results in 



D n (r) 



2™r( 



K+1 N 

2 / 



er„ = const. 



(6) 



where T (x) is the Gamma function, which are equal to 
the central absolute moments of a Gaussian distribution 
N (0, 2cr 2 ) , because the components of u t are (by con- 
struction) Gaussian distributed. 

Our basic idea for introducing spatial correlation be- 
tween the fluid particles is to formulate an Heisenberg-like 
Hamiltonian which is then minimized. In more detail let 
us consider a set of n p fluid particles at positions Xj in 
a cubic volume of linear size L v with periodic bound- 
ary conditions. Every fluid particle i has a velocity u^, 
acceleration and (local) energy dissipation £j. These 
particles are independently integrated in time according 
to Eqs. (I5j). We then define an Hamiltonian 



u = -j2j( r ij) 



u, ■ u, 



(7) 



where ry = |xj — X{| is the distance between particle 
i and j, and J (rij) is a distance dependent coupling 
function. The choice of this coupling function is not known 
a priori and heavily influences the resulting correlations. 
The minimization of the Hamiltonian H is achieved by 
a standard Metropolis algorithm [2S]: A pair of fluid 
particles is randomly selected and a new "configuration" 
is proposed consisting of interchanging the two particles. 
Then the change in energy AE is calculated and the 
interchange is accepted with probability 



p (AE) = mini 1, exp(-AE/T) 



(8) 



where T is an artificial "temperature" , chosen to be "low" . 
This kind of correlation steps basically just aligns the fluid 
particles to minimize the Hamiltonian |7|. This ensures 
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TABLE I. Parameters used throughout the simulations in this 



10 1 



paper. 


Rx 


372 


740 


1115 


crl 

V 


0.4802 m/s 
25.0 m 2 /s 3 
0.000001 m 2 /s 


0.9604 m/s 
25.0 m 2 /s 3 
0.000001 m 2 /s 


1.4406 m/s 
25.0 m 2 /s 3 
0.000001 m 2 /s 


Jo 

a 

r c 


0.00001 

-1.61 

0.01m 


0.00001 
-2.35 
0.02 m 


0.00005 
-3.31 
0.03 m 


L v 
n p 
T 


1.0m 
27000 
0.01 


1.0 m 
27000 
0.01 


1.0m 
27000 
0.01 



that the total kinetic energy and momentum of the system 
are conserved, and additionally also the temporal statistics 
of the SDEs are (on average) unchanged. The correlation 
steps are then included in the whole algorithm, such that 
between two time integration steps of all the SDEs a 
certain number of correlations steps is performed. The 
number of correlation steps is chosen such that after one 
Kolmogorov time r ?) = (^/(e)) 1 / 2 about n p correlation 
steps are performed. 

The coupling function J (r^ ) in Eq. ([7| has to be ad- 
justed such that the desired structure functions are ob- 
tained. Here J (fy) has been chosen as a power law 



J (nj) 



(9) 



Other choices for J (ry) were used as well, but we found a 
power law to give the best results and to be the easiest to 
control. Changing the parameter a leads to a horizontal 
shift of the resulting structure functions and variations 
of Jo change the exponents of D n (r). Since the model 
Hamiltonian ([7]) is purely empirical, we have no concise 
explanation why a power law for J(rij) gives the best 
result. It could be possible that other functional forms 
for the coupling function may give similar results. 

We performed simulations with n p = 27000 fluid parti- 
cles for three different Reynolds numbers at the Taylor 
microscale R\ = 372, 740 and 1115. The parameters used 
in these simulations are listed in Tab. [i] The parameter 
r c is a cut-off distance for the power law (J9j) . One reason 
for such a cut-off is to reduce the complexity of the sim- 
ulation. Calculating the energy in the Hamiltonian Q 
would grow quadratically with the number of particles if 
J (rij) had an infinite range. Introducing a cut-off reduces 
the complexity of this calculation drastically, since only 
particle pairs within a distance r c need to be considered. 
The second reason for using a cut-off is that the infinite 
range of the power law gives rise to too strong correlations 
between the fluid particles at large distances. We found 
that choosing a cut-off distance r c , which is a bit smaller 
than the turbulence length scale L to give the best results. 
Since L increases with increasing R\ , the values of r c are 
different for the two Reynolds numbers. 
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FIG. 1. Measured second order spatial velocity structure 
function D2 (r) for Reynolds numbers R\ = 372, 740 and 
1115. For intermediate distances the desired r 2//3 scaling is 
observed. For smaller distances the structure function does 
not show the behaviour given in Eq. Q, because our model 
is not able to correlate the velocities sufficiently at small 
distances. The insets show the higher order velocity structure 
functions rescaled with the large distance limits from Eq. |6| 
for n = 2, . . . , 8 (top to bottom) for R\ = 740. The increase 
of the exponents of the power laws is clearly visible. 



Fig. [T] shows the second order structure function 
D2 (r /t)) for the three considered Reynolds numbers 
R\ = 372, 740 and 1115. Depending on this Reynolds 
number, we can observe the desired dependence for large 
distances given by Eq. Q. Unfortunately for small dis- 
tances the measured structure functions do not agree with 
the desired functions Q. To understand this behaviour 
we first note that the second order structure function can 
be rewritten [9] as 



where 



D 2 (r) = 2al - 2i? y (r) , (10) 



i?ll(r) = -(u(x + r,i).u(x,t)) (11) 



is the (longitudinal) spatial velocity correlation function. 
The factor 1/3 is introduced because we are only calcu- 
lating the correlation of one component. At this point 
it is important to remark that the Hamiltonian ([7]) only 
correlates the directions of the fluid particles, but not 
their magnitudes. This can be seen by writing 

(r) = R & (r) • R d (r) (12) 

with the mean of the product of the velocity magnitudes 



i4(r)=*(|u(x + r)||u(x)|), 



(13) 



and the spatial correlation function of the directions of 
the fluid particle velocities 



Rd(r) 



u (x + r) • u (x) 
|u(x + r)||u(x)| 



(14) 
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FIG. 2. Measured spatial velocity correlation functions at 
R\ = 740. The velocity correlation 7?|| (r/rf) is too small at 
small distances. First, this is due to the facts that our model is 
not able to perfectly correlate the directions of the fluid particle 
velocities for small r as shown by Rd{r/r/). Second, the 
Hamiltonian |7f does not correlate the magnitudes of the fluid 
velocities which leads to R a {r/rf) being almost independent 
of r. 

Fig. [2] shows measurements of these three functions for 
R\ = 740. From this plot it is clearly visible that the 
correlation function R(r/rf) j<j\ does not approach unity 
in the limit of distance tending to zero, and because of 
that the second order structure function does not go to 
zero as desired. This is a result of two factors: First the 
directions of the fluid particles are not correlated strongly 
enough on small distances, as can be seen also in Fig. [2j 
since R& {r/rj) does not go to one. This behaviour is 
reasonable since there are only a rather small amount 
of fluid particles present in the system and it may not 
be possible to find "a perfect match" for every pair of 
fluid particles. One should note that only a small amount 
of particles (less than 1%) are located close enough to 
even give a contribution to the correlations on these small 
distances. The second reason why the correlation function 
R (r/ri) ja\ does not go to unity is that the magnitudes 
of the velocity vectors are not correlated. This can be 
seen by the fact that i? a (r / rj) is almost independent on 
r. The value of i? a (r/rj) can be calculated analytically: 
The components of the velocity vector of fluid particles 
are basically normally-distributed random variables. As- 
suming for the time being a velocity variance o\ = 1 m 2 /s 2 
results that the absolute values of such vectors follow a 
^-distribution with 3 degrees of freedom and their proba- 
bility density function (PDF) is 

* fe3) =75F(r v * ,/2 - (i5) 

In general, the PDF of the product of two statistically 
independent random variables x and y, with PDF px(x) 
and py (y) can be calculated by 



Pz=x y (z) = 




FIG. 3. Numerically measured exponents of the spatial velocity 
structure functions up to eighth order. For n > 3 the measured 
exponents are smaller than the prediction of Kolmogorov's 
K41 theory. This result is confirmed by experiments [20] , 

For two x-distributed random variables this integral gives 
Pz (z) = 2 z 2 K Q (z) , (17) 

7T 

where Kq (z) is a modified Bessel function of second kind. 
The expectation value of z is then calculated to be 

(z) = -. (18) 

7T 

Only considering the component parallel to r and intro- 
ducing an arbitrary velocity variance a\ gives 

Re.{rh) = ~al. (19) 

We further calculated the higher order velocity struc- 
ture functions D n (r/rj) up to eighth order. For R\ = 740 
these curves are shown in the insets of Fig. [T] For the 
other two Reynolds numbers the curves are similar. At 
large distances the fluid particles are not correlated any- 
more and the values of D n (r/rj) should approach the 
constants given in Eq. The resulting curves are there- 
fore rescaled with these constants. The increase of the 
exponents of the power laws for intermediate distances 
is clearly visible. For larger values of n the power laws 
are less pronounced and it becomes more difficult to de- 
termine the corresponding exponents. We also note, that 
with increasing Reynolds number the transition region 
from the power law to the large scale behavior becomes 
wider and therefore the power law region also becomes 
less pronounced. We measured the exponents £ n of the 
power laws and the results are shown in Fig. |3j together 
with the theoretical value of Kolmogorov's K41 theory 
[H] Cn = f ■ It can be seen that the resulting exponents 
are different for the three considered Reynolds numbers. 
This difference may be reduced by further fine tuning of 
the parameters Jo, a, or r c . The measured exponents 
also move away from Kolmogorov's prediction for higher 
orders — a property well known from experiments [§1 120] . 
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FIG. 4. (Color online.) Measured mean squared displacement 
of fluid particles at R\ = 372, 740 and 1115 rescaled by the 
short time behaviour 3a 2 t 2 . Simulations without correlations 
(solid lines) show a clear ballistic (~ t 2 ) regime at short times 
and a transition to a normal diffusion (~ t) regime for long 
times. The latter is indicated by the solid black line. For higher 
Reynolds numbers this transition happens later since TL is 
larger for higher Rx- For the cases with correlations (dashed 
lines) there is no clear ballistic range visible for R\ = 372 and 
740 and the transition to a diffusion regime happens much 
earlier than in the case without correlations. For R\ = 1115 
the results are closer to the uncorrelated ones. 



III. DISPERSION OF FLUID PARTICLES 

With the model described in the last section we are 
able to generate a turbulent fluid velocity field that (at 
least on not too small scales) obeys the desired spatial 
statistics. We now take a closer look at the dispersion of 
individual fluid particles. 

We start by looking at the mean squared displacement 
(MSD) of single particles. For times much smaller than 
Tl one expects to observe a ballistic scaling and for large 
t normal diffusion should be observed. This means that 
the MSD should behave as 



(|x(i)-x | 2 ) = 



3a 2 t 2 t < T L 
2a 2 u nt t>T L 



(20) 



where Xo is the position of the particle at time t = 0, and 
the brackets (•) here mean averaging over many particle 
tracks. 

In Fig. [4] we show the MSD rescaled by the short time 
behaviour for the three simulated Reynolds numbers with 
and without spatial correlations. For the simulations with- 
out correlations (corresponding to solid lines in the plot) 
clear ballistic scaling for small t is visible. For large t the 
transition to the normal diffusion regime (whose slope is 
indicated by the solid black line) is clearly observed, and 
since Tl is larger for higher R\ this transition happens 
later for the higher Reynolds numbers. For the simula- 
tions with spatial correlations (corresponding to dashed 
lines in the plot) the observed behaviour is different for 
the three Reynolds numbers. For R\ = 372 and 740 



no clear ballistic range is observed and the transition to 
normal diffusion is much faster than in the case without 
correlations. For i?> = 1115 on the other hand the influ- 
ence of the spatial correlations is much less pronounced. 
For short times the t 2 scaling can still be observed and 
the transition to the diffusive regime happens close to the 
desired point at time T^ 115 . 

Due to the interchanging of fluid particles during 
correlation steps the velocity of a certain particle may 
change rather abruptly. This instantaneous modification 
of the fluid particle velocity results, on average, in a 
reduced Lagrangian correlation time Tl, and therefore 
the transition to the normal diffusion regime happens 
earlier in the cases with spatial correlations. Because the 
Lagrangian correlation time T^ 115 > T^ 40 > T 312 already 
established spatial correlations should be preserved longer 
for higher Reynolds numbers. Therefore it is reasonable 
to assume that in a certain time interval more of the 
proposed fluid particle exchanges are accepted for the 
smaller Reynolds numbers, which again implies that the 
Lagrangian correlation time is reduced more for smaller 
Reynolds numbers than for higher ones compared to the 
corresponding correlation times without correlations. 
This may be a reason why the curve of the MSD is less 
influenced by the correlation steps for R\ — 1115. 

Next we investigate the dispersion of pairs of fluid 
particles by measuring their mean squared separation 
(MSS). Let A (t) denote the distance between two fluid 
particles at time t, and Ao the distance of these particles 
at time t — 0. For this relative dispersion one expects to 
observe different scaling subranges as 



^ 2 ((e)A ) 2/3 i 2 ^«t«t 

|A(/)- A„| J ; = {g(e)t 3 t « * « T E 

A<J 2 u T L t i»T L 

(21) 

where to = (Aq/ (e)) 1 ^ 3 and the large eddy lifetime Te = 
a 2 /(e). The first range, which is scaling with t 2 , goes 
back to Batchelor [3J. In this regime particles remember 
their initial velocity differences and therefore their motion 
is highly influenced by the form of the velocity structure 
functions D n (r) . The second regime is known as the 
Richardson t 3 law, even so it was formally first introduced 
by Obukhov [2] 15 years after Richardson's famous work 
PQ about turbulent diffusivity. The constant g is known as 
the Richardson constant, whose exact value is still under 
debate [5). The third scaling range is twice the MSD 
result and reflects the fact, that after a long time t ^> Tl 
the motion of fluid particles is not correlated anymore 
and both particles move independently from each other. 

Fig. [5] shows the measured relative dispersion of fluid 
particle pairs. As reported previously 6, 15 the behaviour 
of the pair separation strongly depends on the value of 
the initial separation Ao- The solid curve corresponds to 
twice the measured MSD of single particles. The vertical 
lines indicate the time scales Te at the corresponding 
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FIG. 5. Measured dispersion of pairs of fluid particles rescaled by t 2 . The expected Richardson t 3 could not be observed. The 
vertical line indicates the large eddy lifetime Te, the arrow shows the range of to in Eq. (21 1, and the transverse lines indicate 
normal diffusion with slope —1. We observe a large dependence on the initial separation between fluid particles Ao for higher 
R\. (a) Results for R\ = 372. Since Te ~ to no Richardson scaling is observed, (b) Results for R\ = 740. For short t a 
Batchelor t 2 scaling is observed. For t in the range of to the increase of the pair separation is slower than t 2 . For the smallest 
initial separation Ao a quadratic scaling can be observed again for times larger then Te. (c) Results for R\ = 1115. Here for the 
smallest initial separation Ao an increase of the pair separation faster than t 2 can be found for Te < t < 1000 t v , which may be 
an indication of Richardson scaling in this region. 



Reynolds number. The values of to range from about 
17 trj at Ao ~ 0.001m up to to ~ 200 t v at large Ao ~ 
0.04 m, which are indicated by the horizontal arrow. The 
transverse lines have slope —1 to identify regions with 
normal diffusion. 

Fig.l5](a) shows the resulting MSS for R\ — 372 rescaled 
with t . Since the eddy lifetime Te is of the same order 
as the values of the time scale to it is n °t surprising, that 
no Richardson t 3 scaling can be observed. Additionally 
only for short times a slight dependence on the initial 
separation Ao can be seen. 

The results for R\ = 740 are shown in Fig. [5] (b) . We 
do not observe Richardson scaling in this case as well and 
therefore the data has been rescaled with t 2 . For short 
times we observe a ballistic regime which is increasing for 
smaller initial separations A . For intermediate times in 
the range of to the pair separation becomes slower and 



tends to a normal diffusion for small initial separation. 
At the smallest initial separations Ao we observed the 
pair separation to almost form a plateau at t ~ 1000 t rp 
which would again correspond to ballistic motion. For 
large times all curves tend to converge to the MSD value 
of single particles. 

In Fig. [5] (c) we show the same measurements for 
R\ = 1115. Again a clear Richardson scaling could not 
be observed. This observation was also reported for ex- 
periments by Bourgoin et al. |15j . The ballistic regime for 
short times is much more pronounced than for R\ = 740. 
For t in the range of Tg-lOOO^ we observe rather differ- 
ent behaviour of the pair dispersion for different values of 
Ao. For Ao > 0.02 m the pair separation is slower than in 
the ballistic case. For Ao ~ 0.01 m we observe a plateau, 
which means that the pair separation is proportional to 
t 2 . For the smallest initial separations Aq < 0.002 m we 
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finally observe an increase in the mean distance between 
pairs of fluid particles which is faster than t 2 . This may 
indicate that a Richardson scaling may be observed in 
this region. 

The systems studied in this paper are quite dilute 
and particles are almost uniformly distributed inside the 
system at any given point in time. We therefore do not 
have any control over the initial separation of particles 
and it is rather difficult to gather enough statistics for 
the MSS at small Aq. Thus the values of Ao were binned 
and the curves of the MSS in Fig. [5] were produced by 
averaging over all pairs within the corresponding bin. The 
binning and subsequent averaging may be a reason for 
not observing a Richardson scaling, because the number 
of pairs with a certain initial distance increases with A 
and therefore the averages may be biased towards the 
upper edge of the bin interval. This bias is particularly 
important for the smallest considered bin where the initial 
separation is most important [Sj for the behaviour of the 
MSS. 



where m, is the mass of particle i, pt the fluid density 
and the drag coefficient Cd is given by 



I 0.44 



(1 + 0.15 Rep 1 



Re p < 1000 
Re p > 1000 



(23) 



This coefficient depends on the particle Reynolds number 
given by 



Re p = — ^ 

v 



(24) 



To quantify the influence of the inertia of particles 
suspended in the fluid usually the Stokes number St is 
used. It is defined by 



St = ^, 

Tf 



where the particle response time r p is given by 



r %Pv 



4 

18 v Pi 



(25) 



(26) 



IV. DISPERSION OF HEAVY PARTICLES 



So far we considered the motion and dispersion of 
fluid particles, which are by definition following the fluid 
streamlines perfectly. We now want to replace the fluid 
particles by heavy real particles and investigate the influ- 
ence of the inertia of these particles on their dispersion. 
Therefore we introduce in addition to the n p fluid par- 



tides the same number of real particles with radius 
density p p and velocity v^. We follow the same idea as 
in Ref. [24] and "glue" one fluid particle to every real 
particle. This means that these two particles are always 
located at the same point in space. Energy dissipation 
£j, fluid particle acceleration a,; and fluid particle veloc- 
ity Uj are still calculated by integrating Eqs. (|5a|H5c), 
but the positions and particle velocities are now 
calculated using the discrete element model (DEM) |27j . 
Particles are treated as soft spheres and particle collisions 
are considered to be elastic with damping. Therefore we 
adopt the linear spring-dashpot model to resolve particle 
collisions. In this model particles are allowed to slightly 
overlap and when they do a repulsive force proportional 
to the overlap is applied. Further explanations of the 
DEM and spring-dashpot model can be found in many 
publications [371 [2H] . The (real) particles are coupled 
to the fluid by an empirical drag law, which gives an 
additional force in the DEM that accelerates the particles. 
For the drag force acting on particle i we use the widely 
used and well established expression 



and Tf is a characteristic fluid time scale. If St is small 
the particles follow the fluid streamlines rather closely, 
and the larger St gets, the less the particle motion is 
influenced by the fluid. In a turbulent velocity field r f can 
range from the small Kolmogorov time t v up to the many 
orders of magnitude larger energy-containing time scale 
Tl- This means that a particle with response time r p is 
only marginally influenced by fluid structures changing 
on time scales Tf < t p . 

In this paper we investigate three different particle 
radii r p = 10~ 5 m, r p = 10~ 4 m, and r p = 10 -3 m. The 
particle and fluid densities are p p = 12.0 kg/m 3 and pf = 
1.2kg/m 3 . This gives particle response times of t p ~ t n , 
t p w 100 t rjl and t p k, 10000 t v . For the smallest particle 
size the collisions between particles were ignored, since 
the particle volume fraction of this system is extremely 
small. In this section we only present results for the two 
higher Reynolds numbers R\ = 740 and 1115. 

Fig. [6] shows the measured MSD of single particles 
rescaled with t 2 and compared to the results for the 
fluid particles. The short time behaviour of Eq. (20) 



contains the mean squared fluid particle velocity 3erJ 
This suggests that in Fig. [6] the curves should go towards 
the mean squared (real) particle velocities Vq (r p ) for 
t — > 0, and indeed we can observe a mean squared velocity 
which is decreasing with r p . This is perfectly reasonable 
since within time intervals of length about Tl , larger and 
therefore heavier particles can be accelerated on average 
only to smaller velocities Vj before the fluid velocity 
changes again. Another observation we can make is that 
the ballistic range for small t is increasing for larger r p , 
which is a result of the increasing particle response time 



jidrag 
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m, 
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(22) 



The smallest particles with r p = 10~ 5 m are of the same 
size as the Kolmogorov length 7/ and therefore one expects 
that the MSD of these particles is almost the same as the 




FIG. 6. Numerically measured mean squared displacement for real heavy particles at two different Reynolds numbers compared 
to the MSD of fluid particles. The main effect of the inertia of the real particles are a reduction of the mean particle velocity 
and a prolongation of the ballistic range at small t. (a) Results for R\ — 740. (b) Results for R\ = 1115. For r p — 10 -5 m the 
MSD of fluid and real particles are almost the same. For r p = 10 -4 m an increase of the ballistic range at short times is visible. 
For long times the MSDs are again almost the same. For the largest particles at R\ = 1115 the ballistic range is increased 
further and the transition to normal diffusion happens later than for the fluid particles. 



one of fluid particles. For larger, and therefore heavier 
particles their response time r p is increased and thus the 
influence of the particle inertia should be visible by an 
enlarged t 2 range of the MSD, because the larger particles 
are less influenced by the fast and small scale velocity 
fluctuations of the fluid velocities and need more time to 
adjust their own velocity to the one of the fluid. Both 
of these effects can be seen in Fig. [6] for R\ — 740 and 
Rx = U15. 

Finally we again measured the dispersion of pairs of 
particles. The resulting MSS for two Reynolds numbers 
R\ = 740 and R\ = 1115 are shown in Fig. [7] As in 
the case of fluid particles we did not observe a clear 
Richardson scaling and therefore the results were rescaled 
by t 2 . The main influence of the inertia of the particle can 
be seen by an increase of the Batchelor t 2 range at short 
times with increasing particle radius. Also the decrease 
of the mean particle velocity vq (r p ) for larger r p can be 
seen by a reduction of the MSS for short times compared 
to the case of fluid particles. For the smallest particles 
with radius r p = 10~ 5 m (Fig. 0(a) and (b)) the resulting 
MSSs are very similar to the one of fluid particles shown 
in Fig. [5] As in the case of fluid particles an increase of the 
MSS faster than t 2 can only be observed for R\ = 1115. 
It is also interesting to note, that the Richardson regime 
at R\ — 1115 is more pronounced for larger particles. 
By looking at the MSD of single particles in Fig. 6] (b) 
one can see that for the largest particles the ballistic 
regime is actually larger than it is for smaller particles. 
This means that the integral timescale Tl is larger in this 
case, and since Tl and Te are usually of the same order 
of magnitude, one can expect that the time interval in 
Eq. (21 1, where a Richardson regime is expected, is larger 



"more time" to show a Richardson i 3 law in the MSD. 



for larger particles than it is for smaller particles. The 
prolongation of this interval then gives heavy particles 



CONCLUSION 



In this paper we presented a new method to model the 
dispersion of particles in turbulence. We used a set of 
SDEs to simulate the temporal evolution of Lagrangian 
tracer particles and introduced spatial correlations be- 
tween them by minimizing an Heisenberg-like Hamilto- 
nian. With this model we were able to produce turbulent 
velocity fields that obey the measured temporal statistics 
and show the correct spatial velocity structure function 
on distances r > 500 77. Investigations of the MSD of 
single particles show a ballistic regime for short times 
and a transition to a normal diffusion regime for large 
times. The algorithm for introducing spatial correlations 
shifts the moment in time where this transition occurs 
toward shorter times and this influence is smaller for 
larger Reynolds numbers. Further the MSS of pairs of 
particles were investigated. We were able to observe a 
Batchelor t 2 range for short times, but failed to see a clear 
Richardson t 3 scaling. Only for higher Reynolds numbers 
and small initial separations indications of a Richardson 
scaling could be noticed. The effects of the inertia of 
heavy particles on their dispersion was investigated as 
well. For short times an increase in the ballistic regime in 
the MSD as well as an increase in the Batchelor regime 
in the MSS could be observed with increasing particle 
radius. Indications of Richardson scaling were observed 
as well. 
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FIG. 7. Numerically measured MSS for different particle radii and Reynolds numbers, (a) & (c) Results for Rx = 740. (b) & (d) 
Results for Rx = 1115. For the smallest particle radius the curves are almost the same as for fluid particles. For larger particles 
the t 2 regime at short times is enlarged and the mean particle velocities are reduced. For the largest particles also the long time 
behaviour of the MSS is influenced. 
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